NONLINEAR AEROELASTIC ANALYSIS USING A TIME-ACCURATE 
NAVIER-STOKES EQUATIONS SOLVER 


Geojoe Kuruvila 1 , Robert E. Bartels 2 , Moeljo S. Hong 3 , and Kumar G. Bhatia 4 

'Flight Sciences Technology, The Boeing Company 
5301 Bolsa Avenue, MC H013-B308, Huntington Beach, CA 92646, USA 
e-mail: geojoe.kuruvila@boeing.com 

2 Aeroelasticity Branch, NASA Langley Research Center 
MS 340, Hampton, VA 23681, USA 
e-mail: r.e.bartels@larc.nasa.gov 

Enabling Technology and Research, The Boeing Company 
PO Box 3707, MC 67-LF, Seattle, WA 98124, USA 
e-mail: moeljo.hong@boeing.com 

4 Aeroelasticity and Multidisciplinary Optimization, The Boeing Company 
PO Box 3707, MC 03-KR, Seattle, WA 98124, USA 
e-mail: kumar.g.bhatia@boeing.com 


Keywords: Nonlinear Aeroelasticity, Limit Cycle Oscillation, Control Surface Freeplay, 
Computational Fluid Dynamics. 

Abstract. A method to simulate limit cycle oscillation (LCO) due to control surface freeplay 
using a modified CFL3D, a time-accurate Navier-Stokes computational fluid dynamics (CFD) 
analysis code with structural modeling capability, is presented. This approach can be used to 
analyze aeroelastic response of aircraft with structural behavior characterized by nonlinearity 
in the force verses displacement curve. A limited validation of the method, using very low 
Mach number experimental data for a three-degrees-of-freedom (pitch/plunge/flap deflection) 
airfoil model with flap freeplay, is also presented. 

1 INTRODUCTION 

Classical aeroelastic analysis methods assume linear aerodynamic and structural behavior. 
While these methods have been quite successful in predicting aeroelastic phenomenon such as 
flutter, their limitations are also well known. They fail to predict instabilities caused by 
nonlinearities in the aeroelastic system. In the past two decades considerable experimental and 
analytical efforts have been devoted to understand nonlinear aeroelasticity. However, accurate 
prediction of instabilities such as LCO, transonic flutter, and buzz remains a challenge [1]. It 
is widely believed that these instabilities are caused by nonlinearities in the fluid, or the 
structure, or both. LCO, for example, is thought to be caused by nonlinear flow due to large 
shock motion or flow separation, or structural nonlinearities arising from freeplay in hinges 
and linkages of control surfaces or material behavior [2], Aging aircraft and combat aircraft 
carrying heavy external stores are known to experience LCO, typically due to structural 
nonlinearities. Aerodynamic nonlinearities can be modeled using Navier-Stokes equations for 
flow simulation. For complex structures like that of an aircraft, it is virtually impossible to 
determine the sources of all structural nonlinearities. Generally a particular source of 
nonlinearity, such as control surface freeplay, is explicitly modeled as part of the fluid- 
structure interaction problem in order to predict system instabilities. 
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About a decade ago, a study of LCO due to control surface freeplay was conducted by 
Conner, Tang, Dowell, and Virgin using a three-degrees-of-freedom aeroelastic wing model 
with piecewise linear flap stiffness [3, 4]. The study included testing in the Duke University 
low-speed wind tunnel and numerical simulations. A system of piecewise linear state-space 
models was used for the numerical simulations. The equations of motion were integrated in 
time using a Runga-Kutta algorithm in conjunction with Henon’s method to locate the 
“switching points” in the flap stiffness as it changed from one state to another. Later, Tang, et 
al. [5] and Gordon [6] successfully used the describing function approach to model this 
nonlinear problem. To model the aero-structure interaction, Conner and Gordon used linear 
potential flow theory based Theodorsen’s method while Tang used Peters’ method. The 
results of their low Mach number (0.02 < M < 0.07) simulations correlated reasonably well 
with experimental data. 

The present study explores the feasibility of using a nonlinear aero-structure interaction code 
for simulating LCO due to control surface freeplay. Such a capability is necessary to 
understand this phenomenon particularly in the transonic regime where the aerodynamics is 
highly nonlinear. This paper describes a method where structural nonlinearity due to control 
surface freeplay is modeled in a nonlinear aerodynamics analysis code. This method has been 
implemented in CFL3D [7, 8], a time-marching Reynold’s Averaged Navier-Stokes code with 
linear structural modeling originally developed by NASA Langley Research Center. As a first 
step, validation of this method has been done using experimental data, from Conner, et al., for 
the three-degrees-of-freedom airfoil model with freeplay. Simulations of this model using 
Navier-Stokes and Euler equations and the comparisons of results with experimental data are 
also presented in this paper. 

2 MODELING CONTROL SURFACE FREEPLAY 

The standard CFL3D simulates aeroelastic phenomenon using the Euler or the Navier Stokes 
equations for modeling aerodynamics and a modal dynamic model for structure [7, 8]. The 
structural model is limited to simulating linear structural behavior. Modeling of structural 
nonlinearity characterized by piecewise linear stiffness variations, commonly seen in control 
surfaces with freeplay, is the focus of this study. An idealized behavior of a control surface 
with such a freeplay is shown schematically in Figure 1. Within a range of deflection angles 
(- S 0 < 8 < S 0 ) the control surface experiences no resistance to its movement. Outside this 
critical range it experiences a linear restoring torque. 



o.o 

Control Surface Deflection 


Figure 1: An idealized depiction of symmetric control surface freeplay 
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The aeroelastic equations of motion for a linear structural model written in vector notation are 
shown in equation (1), where m , d , and k are the mass, damping, and stiffness matrices, 8 
the physical displacements vector and f the aerodynamic forces vector. 

m8 + d8 + kS = f (1) 

When some modes have freeplay like the one shown in Figure 1, equation (1) can be modified 
as, 


m5 + dS + k5 + k8 = f (2) 

where, for modes with freeplay, i , 
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and, for modes without freeplay, j , £ . = 0 . Note that equation 3 can be easily modified to 
model asymmetric freeplay behavior. 
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Figure 2: Decomposition of freeplay behavior (a) into linear (b) and residual (c) components 


In equation (2), the restoring torque due to freeplay is modeled as the sum of a linear term 
(k5) and a nonlinear residual term (ks). This decomposition is shown graphically in Figure 2. 
After moving the residual term to the right hand side and performing modal transformation, 
equation (2) in generalized form is, 

Mq + Dq + Kq = F-(|) T k8 (4) 


where, M , D , K , q , and F are the generalized mass, damping, stiffness, displacements, and 
aerodynamic forces, respectively. They can be expressed as follows. 

M = (|) T m(|) 

D = <|> T d<t) 

K = <|) T k<t) (5) 

8 = <(>q 
F = <|> T f 
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where <j> are the eigenvectors of the system . mS + kS = 0 . <|) T kS are the generalized residual 
forces. In the absence of freeplay, these terms would be zero. 

In CFL3D the modal equations of motion are cast in the state-space form and integrated in 
time. For simplicity, the off-diagonal terms in the damping matrix are assumed to be zero. It 
results in a 2x2 block diagonal system matrix that is easily integrated using the Laplace 
transform technique. For modeling freeplay, the residual terms due to nonlinearity in stiffness 
are added to the right-hand-side of the state-space equations. The left-hand-side remains the 
same and therefore the technique for integrating the state space equations remains unchanged. 

3 METHOD VALIDATION MODEL 

Data from experiments conducted by Conner are used to validate the formulation described 
above. A schematic of the model used is shown in Figure 3. It is a two-dimensional 
rectangular wing with a full-span trailing-edge flap. The wing airfoil cross-sections are NACA 
0012 with the aft quarter of the airfoil forming the flap with the hinge at % chord location. 
The model has three-degrees-of-freedom, namely, pitch, plunge, and flap deflection. The pitch 
axis is at !4 chord location. The setup allowed the freeplay band to be an adjustable parameter. 
For completeness the parameters of the experimental model are given in Table 1. 



The displacement vector for this model can be written as, 

8 = {h a pf (6) 

where h , a , and [3 are plunge, pitch, and flap deflection, respectively. The flap can be preset 
to have a desired freeplay zone of ±/? 0 , i.e., between - /? 0 < f3 < /? 0 its stiffness is zero. The 
generalized residual term on the right hand side of equation (4) becomes 

<t> T kS = <|) T kjo 0 pf (7) 


where /? is given by equation (3). 
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Geometric parameters 

Chord 

0.254 m 

Span 

0.52 m 

Semi-chord, b 

0.127 m 

Elastic axis, a w/r/t b 

-0.5 

Hinge line, c w/r/t b 

0.5 

Mass (inertial) parameters 

Mass of wing 

0.62868 kg 

Mass of aileron 

0.18597 kg 

Mass/length of wing-aileron 

1.558 kg/m 

Mass of support blocks 

0.47485 x 2 kg 

S a (P er span) 

0.08587 kg m 

Sp (per span) 

0.00395 kg m 

I a (per span) 

0.01347 kg m 2 

Ip (per span) 

0.0003264 kg m : 

Stiffness parameters 

K a (per span) 

37.3 kg m/s 2 

K p (per span) 

3.9 kg m/s 2 

K h (per span) 

2818.4 kg/m/s 2 

Damping parameters 

Ga (log-dec) 

0.01626 

(log-dec) 

0.0115 

Ch (log-dec) 

0.0113 


S a , Sp\ Static moments of wing per unit length of wing-aileron 
(relative to a) and aileron (relative to c ), respectively. 

I a , Ip: Moment of inertia per unit length of wing-aileron and 

aileron about a and c, respectively. 

Table 1: Parameters of the experimental model [3, 4, 5] 

4 RESULTS 

4.1 Characteristics of the linear system 

From the mass, inertial, stiffness, and damping parameters of the experimental setup, given in 
Table 1, the coupled natural frequency of each degree-of-freedom is determined. Comparison 
of these results with those of Conner’s is shown in Table 2. The numerical results agree with 
one another. The plunge and pitch frequencies agree reasonably well with experiment. The 
flap frequency shows a 4% discrepancy with the experiment. Table 3 shows a similar 
comparison of linear flutter speed and frequency. Here the discrepancies with experiment are 
larger. CFL3D result shows the least difference. Reference 3 mentions degradation of the 
stiffness of the springs, especially the pitch-spring, due to fatigue. This could be a reason for 
the discrepancies. 
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Coupled Natural 
Frequency (Hz) 

Experiment 

Current 

Numerical 

Conner’s 

Numerical 

Plunge 

4.375 

4.448(1.7) 

4.455 (1.8) 

Pitch 

9.125 

9.197 (0.8) 

9.218(1.0) 

Flap Deflection 

18.625 

19.387 (4.1) 

19.442 (4.4) 


Table 2: Natural frequencies of the linear system (% error relative to experiment in parenthesis) 


Linear Flutter 

Experiment 

CFL3D 

Navier-Stokes 

Conner’s 

Numerical 

Speed (m/s) 

20.6 

22.6 (9.9) 

23.9(16.0) 

Frequency (Hz) 

5.47 

5.87 (7.3) 

6.112(11.7) 


Table 3: Flutter speed and frequency of the linear system (% error relative to experiment in parenthesis) 


4.2 Navier-Stokes solutions with freeplay using measured stiffness 

Results of Navier-Stokes analyses with measured stiffness, for freeplay of ±2.12°, are shown 
here. A three-dimensional grid, but with only one cell across the span, is used for these 
simulations. The grid has 353x97 points in the streamwise and normal directions, respectively, 
with 255 points on the airfoil. The freestream velocity for these analyses varies between 7 and 
23 meters per second. Sea-level conditions are assumed to determine the other flow 
parameters. The Reynolds number based on chord is 520,000. Note that at this Reynolds 
number, only about 60% of the airfoil is likely to have turbulent flow. However, fully- 
turbulent Navier-Stokes simulations using the Spalart-Allmaras turbulence model are 
conducted. Each analysis is run for 200,000 time steps with 1 0 sub-iterations. One time-step is 
equal to 0.1 millisecond, making each simulation 20 seconds long. 

The overall comparisons of Navier-Stokes simulations with experimental results are shown in 
Figure 4. The scale, U /U f referred to as the speed-ratio, is the ratio of the flow speed to the 

linear flutter speed. During the experiment, as the speed was increased, five types of system 
behavior were observed. At first the initial perturbations damped out. Then, between speed- 
ratio of 0.23 and 0.37, the system displayed a simple low frequency behavior. From 0.37 to 
0.48, the system lost periodicity. Between 0.48 and 0.55 the system became periodic again 
and the flap oscillated in a “complex” low frequency limit cycle. Beyond speed-ratio of 0.55, 
the system was abruptly attracted to a high frequency limit cycle. The speed ranges of these 
five types of behavior are indicated by the lower color bars in Figure 4. 

The upper color bars indicate the system behavior simulated by CFF3D. Overall the stages 
described by the experiment are also seen in the CFD simulation. However, there is 
considerable discrepancy in speed ranges at which they occur. In these simulations the simple 
oscillations do not appear until the speed-ratio reaches 0.39. When it reaches 0.42, the 
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Figure 4: Comparison of Navier-Stokes simulations with experiment p 0 = ±2.12° 


oscillations lose periodicity and remains in this state until it becomes 0.53 when the 
oscillations become periodic again and the flap moves in a “complex” low frequency limit 
cycle. At a speed-ratio of about 0.58 the system becomes non-periodic again until it reaches 
0.63 when it begins to exhibit high frequency LCO. The non-periodic behavior, indicated by 
the orange bar, is not specifically mentioned in the description of the system behavior during 
the experiment. 

The relative power spectral density for flap frequencies up to 20 Hz at different speed-ratios, 
shown in Figure 5, presents another view of the system behavior. At first the initial 
perturbations damp out (purple). Then, as the speed increases (0.39-42), frequencies between 
4 and 5 Hz begin to dominate (blue). This results in the “simple periodic low frequency 
LCO.” At higher speeds (0.42-0.53), in addition to the 5 Hz, 10 Hz oscillations also begin to 
appear (green). It is weaker than the 5 Hz one, but gains strength as the speed increases. 
Together they produce a mildly non-periodic behavior. Between speed-ratios of 0.53 and 0.58 
(yellow), the strengths of these frequencies stabilize and the oscillations become periodic 
again. The flap motion has the “complex” structure due to the presence of two dominant 
frequencies. Beyond 0.58 (orange), the 5 Hz mode loses strength rapidly while the 10 Hz 
mode gains strength. When the speed-ratio reaches 0.63 (red), the 5 Hz mode dies completely 
and only the 10 Hz mode remains resulting in the high frequency limit cycle. Beyond this, 
only one dominant mode exists and its frequency increases with speed. It reaches about 12 Hz 
near the linear-flutter speed. The power spectral density of plunge and pitch modes are 
relatively small compared to that of the flap oscillation mode. 
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Figure 5: Relative power spectral density of flap oscillations for J3 0 = ±2.12° 


For quantitative comparison of the CFL3D results with Conner’s results, time histories of 
displacements of the three modes are normalized as follows 



where h, a, and p are the time histories of plunge, pitch, and flap displacements, 
respectively. H , A, and B are the corresponding non-dimensional displacements and b is the 
semi-chord. From the non-dimensional amplitudes and using the time step At, the r.m.s. 
values are calculated as, 


Am Pr.m.s 



( 9 ) 


Comparisons of CFL3D Navier-Stokes solutions with Conner’s results are shown in Figures 6 
and 7. His results include experimental data as well as results from numerical simulation. The 
aeroelastic model used in his numerical simulations is the linear potential equation based 
Theodorsen’s method. Two linear models based on this method are combined to capture the 
nonlinear (piecewise linear) nature of the freeplay. The numerical model appears to scale with 
the size of freeplay. Experimental data for three freeplay bands ( Po = ±2.12°, ±1.83°, and 
±1.15°) are presented. Note that there is significant scatter in the experimental data, especially 
at low speed-ratios where the LCO has multiple dominant frequencies. 
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Figure 6: Comparison of normalized r.m.s. amplitudes of plunge and pitch modes 



Figure 7: Comparison of normalized r.m.s. amplitude of flap deflection and dominant flap frequency 


Figures 6 and 7 show that, Conner’s results and CFD simulations have similar trends. 
However, their magnitudes have considerable discrepancies. Normalized r.m.s. amplitude of 
plunge at lower speed-ratios are closer to Conner’s numerical results than the experimental 
data. At higher speeds, there is generally better agreement with Conner’s results. For the pitch 
mode, CFL3D predicts the general trend but the discrepancy at higher speeds is large. 
Normalized r.m.s amplitude of the flap mode shows reasonable agreement, but at higher 
speeds the difference is large. The overall agreement of the flap LCO frequencies is 
reasonable. As mentioned before, at very low speed-ratios the CFL3D results do not show any 
LCO, while Connor’s experimental and numerical results show simple periodic LCO. The 
reasons for the discrepancies are not entirely clear. Unknown inconsistencies between 
experimental and CFD parameters could be one. Based on other simulations (presented 
below), there is indication that dissipation in flow modeling may also have a role in causing 
the deviation. 

4.3 Euler solutions with freeplay using measured and adjusted stiffness 

Results of simulation conducted to evaluate the sensitivity to fluid viscosity and structural 
stiffness are presented. In one set of numerical experiments the cases presented in section 4.2 
are repeated using Euler analysis. In the next set, the cases are repeated with Euler analysis 
and adjusted stiffnesses. The stiffnesses are adjusted to match the measured coupled natural 
frequencies of the experimental setup. The presumption here is that the frequencies, in 
general, can be measured more accurately than stiffnesses. The values of measured and 
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adjusted reduced stiffnesses are shown in Table 4. Figures 8, 9 and 10 show the results of 
these analyses. 


Stiffness 

Measured 

Adjusted 

Plunge (kg/m/s 2 ) 

2818.4 

2252.9 

Pitch (kg m/s 2 ) 

37.3 

37.1 

Flap (kg m/s 2 ) 

3.9 

3.5 


Table 4: Measured and adjusted stiffness 




Figure 8: Comparison of normalized r.m.s. amplitudes of plunge and pitch modes 




Figure 9: Comparison of normalized r.m.s. amplitude of flap deflection and dominant flap frequency 

Figure 8 show that Euler analysis moves the results closer to the inviscid results of Connor’s. 
For the plunge mode, this shift further increases the discrepancy between experimental and 
CFD results. For the pitch mode, at higher speeds, the Euler results falls within the scatter of 
experimental data. For the flap mode (Figure 9), there are no major changes at the lower 
speeds. At higher speeds, the r.m.s. of the amplitudes are slightly lower than Navier-Stokes 
solutions. For plunge and pitch mode, the changes in results when using adjusted verses 
measured stiffness are less than the changes from Navier-Stokes to Euler solutions. For the 
flap mode, at the higher speeds, the changes in amplitude are relatively larger. The differences 
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in solutions (not shown) between Navier-Stokes solutions with measured and adjusted 
stiffness values are analogous to the difference between similar Euler solutions. 

A summery of results of Navier-Stokes and Euler analyses, using measured and adjusted 
stiffnesses and their comparison with experiment are shown as histograms in Figure 10. A 
noticeable change when using Euler analysis is that LCO appears at lower speeds than in the 
case of Navier-Stokes solutions. Lower dissipation in the Euler formulation may be helping to 
capture LCO at the lower speeds. The regions with “simple periodic low frequency” LCO and 
“high frequency” LCO are wider than Navier Stokes results and the other regions are 
narrower. Adjusting stiffnesses to match the measured frequencies appear to marginally move 
the numerical results closer to the experimental data. Overall, it is not clear if switching to 
Euler solutions and/or using adjusted stiffness improved the correlation with experimental 
data. 
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Figure 10: Comparison between experiment, Navier-Stokes and Euler results 


4.4 Sensitivity to width of freeplay band 

Simulations with a freeplay band of ±1.15° are compared to those with freeplay of ±2.12° in 
Figures 11 and 12. Euler analyses with stiffness adjusted to match experimentally measured 




Figure 1 1 : Sensitivity of plunge and pitch amplitudes to width of freeplay region 
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Figure 12: Sensitivity of flap deflection and dominant flap frequency to width of freeplay region 

frequencies are used here. While Conner’s numerical results scale with the freeplay band, as 
defined by equation (8), CFL3D results do not. Also note that the experimental data do not 
scale well either. 


4.5 Sensitivity to time-step size 

Results of a time-step sensitivity study are shown in Figures 13 and 14. All the solutions 
presented earlier use a constant time-step of 0.1 millisecond. The flap deflection mode of the 
experimental setup has the highest coupled natural frequency and it is about 20Hz. In the CFD 
simulation, one cycle of the highest frequency is discretized into about 500 time intervals. In 
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Figure 14: Variation of amplitudes and frequency with varying time-step size, U/Uj=0.5 


this time-step study, selected cases are repeated with time-steps of 0.2, 0.5 and 1.0 
millisecond, i.e., each cycle of the highest frequency is modeled with 250, 100, and 50 points. 
Figure 13 shows the responses of the flap mode undergoing “complex” periodic low 
frequency LCO. It can be seen that the amplitude and the phase are sensitive to time-step size. 
Figure 14 shows the changes in amplitudes of the three modes and that of the flap frequency. 
It also shows the percent change with respect to the solution with the smallest time-step size. 
As the time-step size increases, the plunge and pitch modes show large difference in 
amplitudes. The change in flap amplitude and frequency are relatively small. Similar 
sensitivities are observed at other flow speeds. A reason for this could be the accumulation of 
error in numerical integration caused by missing the “switching point” of the piecewise linear 
freeplay behavior. This is identified as a key issue by Conner who uses an adaptive time- 
stepping algorithm to accurately integrate to the comer of the two linear regions. However, 
from a CFD time-accuracy point of view, it is often necessary to use a time-step smaller than 
0.1 millisecond, especially for complex configurations. So, special treatment to capture the 
switching-point may not be worthwhile. 

5 CONCLUSIONS 

A method to simulate LCO due to control surface freeplay has been formulated and 
implemented in CFL3D, a time-accurate Navier-Stokes code. This method can be used for 
simulation of LCO of complete aircraft that has freeplay in one or more control surfaces. A 
reasonable validation of the method has been accomplished using very low Mach number 
experimental data for a simple three-degrees-of freedom system. There is indication that 
excessive dissipation in the CFD method may be inhibiting LCO at the low end of the 
experimental speed range. Uncertainties in the measured stiffnesses and possibly the damping 
coefficients (not investigated here) may also be contributors to the discrepancies between 
simulation and experiment. 

An adaptive time-stepping scheme to integrate accurately to the “switching point” is not 
helpful here as very small time-steps are needed for time-accuracy of the aerodynamic 
solution. While the small time-step requirement is not a major constraint for this 2- 
dimensional problem, it can be an impediment to analyzing complete aircraft configurations. 
Finally, the grid perturbation scheme in CFL3D worked flawlessly for this geometrically 
simple problem. However for complex 3-dimensional configurations, robust grid perturbation 
remains a challenge. Further analyses using aircraft configuration at transonic Mach numbers 
need to be performed to fully evaluate this approach for use in industrial applications. 
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